R Notebook

\(V_{kn}\) functions and Prior cluster distribtion.

Data generation

Gaussian mixture:

\[Y \sim \frac{1}{2}N(1,0.2) +\frac{1}{2}N(5,0.2)\]

Model

\[\begin{equation} \begin{aligned} (X_i \mid (\mu_i,\tau_i)) &\stackrel{iid}{\sim} N(\mu_i,\tau_i^{-1})\\ ((\mu_i, \tau_i) \mid G) &\stackrel{iid}{\sim} G \\ G &\sim \mathcal{P} . \end{aligned} \end{equation}\]

Firstly, assume \(\tau_i=1\),then we have \(X_i \sim N(\mu_i, 1)\), base measure \(G_0 = N(0,1)\), \(Y_i= (\mu_i)\). For sampling we use Polya Urn characterization

\[\begin{equation} p(Y_i | {Y_{-i}}) = \frac{\alpha}{\alpha+n-1} G_0 + \frac{1}{\alpha + n -1} \sum_{j=1}^{n-1} \delta_{Y_j} \qquad \forall i=1,\ldots ,n. \end{equation}\] Gibbs sampler: \[\begin{equation} \begin{aligned} G &\mid \theta, X \sim DP(\alpha +n,\frac{\alpha}{(\alpha + n)} G_0 + \frac{1}{\alpha + n)}\sum_j \delta_{\theta_j)}\\ \theta_i &\mid \theta_{-i},X, G \sim \sum_{j\neq i}q_{ij} \delta_{\theta_j} + r_i G \\ q_{ij} &= b N(x_,\theta_j)\\ r_i &= b\alpha \int f(x_i, \theta)dG_0(\theta), \text{Here: } N(x_i, \theta,1)N(\theta, 0,1)\propto N(x_i/2, 1/2).\\ \end{aligned} \end{equation}\]

Prior

Gbase<- function(num, mean=0, sigma=1){
  return( rnorm(n=1, mean, sigma))
}
PUrn_sample<-  function(n, alpha, Gbase){
  #theta.1 sample from base measure
   theta<- Gbase(1)
   group_vec <- c(1)
   for (i in 2:n){
     #probability to sample from base measure
     p1<- alpha/ (alpha + i-1)
     #probability to sample from existing atoms
     p2<- 1/ (alpha + i-1)
     pvec<- c(p1, p2*group_vec)
    # prob <- c(alpha/(length(count.vector)+alpha), count.vector/(length(count.vector)+M))
    group_num<-  sample(0:length(group_vec), prob=pvec, size=1)
    if (group_num == 0){
      #add new atom
      theta <- c(theta, Gbase(1))
      group_vec <- c(group_vec, 1)
    }else{
      #update existing group for group_vec
      group_vec[group_num] = group_vec[group_num] + 1
    }
  }
  return(list(theta=theta, group= group_vec))
}
alpha_0<-1
n<- 500
x<- PUrn_sample(n, alpha_0,Gbase)
k_groups<- x$group
weights<- c(k_groups/(alpha_0+n), alpha_0/(alpha_0+n))
  
df_weights <- data.frame(matrix(NA, nrow =length(weights), ncol =1))
df_weights$pw<-weights
df_weights$tr<-1:length(weights)
pl_weigths<- ggplot(df_weights, aes(x=tr, y=pw)) +
  geom_segment( aes(x=tr,xend=tr,y=0,yend=pw)) +
  geom_point( size=0.5, color="red", fill=alpha("blue", 0.3), alpha=0.4, shape=21, stroke=2)+  labs(title=paste0("Prior weights"))+
  theme_bw() + theme(axis.text.x = element_text(angle = 0, hjust = 1,size = 10), strip.text = element_text(size = 15),legend.position = "top", plot.title = element_text(hjust = 0.5))
pl_weigths
df_weights <- data.frame(matrix(NA, nrow =length(weights), ncol =1))
df_weights$pw<-weights
df_weights$tr<-c(x$theta,20)
pl_weigths<- ggplot(df_weights, aes(x=tr, y=pw)) +
  geom_segment( aes(x=tr,xend=tr,y=0,yend=pw)) +
  geom_point( size=0.5, color="red", fill=alpha("blue", 0.3), alpha=0.4, shape=21, stroke=2)+  labs(title=paste0("Prior weights"))+
  theme_bw() + theme(axis.text.x = element_text(angle = 0, hjust = 1,size = 10), strip.text = element_text(size = 15),legend.position = "top", plot.title = element_text(hjust = 0.5))
pl_weigths

 PolyaUrnGenerator <- function(num, M, G0.generator) {
   X.vector <- G0.generator(n = 1)
   count.vector <- 1
   for (i in 2:num){
     prob <- c(M/(length(count.vector)+M), count.vector/(length(count.vector)+M))
     category <- which(rmultinom(n=1, size=1,prob=prob)==1) - 1
     if (category == 0){
       X.vector <- c(X.vector, G0.generator(n = 1))
       count.vector <- c(count.vector, 1)
     }else{
       count.vector[category] = count.vector[category] + 1
     }
   }
   return(list(X.vector, count.vector))
 }
## Define a Polya Urn sampler
PolyaUrnSampler <- function(n, M, cutoff, num = 10000, G0.generator=rnorm){
  output <- matrix(NA, nrow = n, ncol = length(cutoff) + 1)
  cutoff <- c(-Inf, cutoff, Inf)
  for(i in 1:n){
    PU.sample <- PolyaUrnGenerator(num=num, M=M, G0.generator=G0.generator)
    for(j in 1:ncol(output)){
      output[i,j] <- sum(PU.sample[[2]][intersect(which(PU.sample[[1]] > cutoff[j]), which(PU.sample[[1]] <= cutoff[j+1]))])/num
    }
  }
  return(output)
}
## 
x<- PolyaUrnGenerator(1,1,rnorm)
y<- PUrn_sample(1, 1, Gbase)

ks.test(f_approx$ecf, ecdf(X))

    One-sample Kolmogorov-Smirnov test

data:  f_approx$ecf
D = 0.79, p-value = 8.882e-16
alternative hypothesis: two-sided

Posterior_P_Urn_exact<- function(X,theta,i,alpha){
  q_weights<- dnorm(X[i], theta)
  r_weight<- alpha*dnorm(X[i], mean=0, sd=sqrt(2))
  #sum should be one
  num_norm<- sum(q_weights)+r_weight
  q_weights_norm<- q_weights/num_norm
  r_weight_norm<- r_weight/num_norm
  label <- sample(0:length(theta),prob=c(r_weight_norm,q_weights_norm),size=1,replace=TRUE)
  if (label==0){
    theta_new<- rnorm(1, mean=X[i], sd=sqrt(1/2))
  }
  else{
    theta_new <- theta[label]
  }
}


sample<- gen_data(mu_vec=c(1,5),sigma_vec=rep(0.2, 2),ns=50,prob_vec=c(0.5,0.5))
X<- sample
theta.matrix <- matrix(NA, nrow = 10000, ncol = length(X))
alpha<-1
#weights.matrix <- matrix(NA, nrow = 100, ncol = length(X)+1)
theta.init<- runif(length(X),-3,3)
theta.matrix[1,] <- theta.init
for (i in 2:nrow(theta.matrix)){
    theta_last<- theta.matrix[i-1,]
    N <- length(X)
    temp_theta <- c(theta_last, rep(NA, N)) 
    for (j in 1:N){
      temp_theta[N+j] <- Posterior_P_Urn_exact(X, temp_theta[(j+1):(j+N-1)],j,alpha)
    }
    theta_S<- temp_theta[(N+1):(N*2)]
    theta.matrix[i,] <-theta_S
}

#burnin period
theta.final <- theta.matrix[5000:10000,] 
 
xgrid <- seq(from= -10, to=2,length=50)
fgrid <- NULL

fbar.H <- function(xgrid,wh)
{ ## return a draw F ~ p(F | ...) (approx)
  
  fx <- rep(0,length(xgrid))
  for(h in 1:length(wh))
    fx <- fx + dnorm(xgrid,m=wh[h],sd=1)
  return(fx/length(wh))
}


gibbs.H <- function(n.iter=500){
xgrid <- seq(from= -10, to=10,length=50)
fgrid <- NULL
plot(density(X),xlab="X",ylab="Y",bty="l",type="l",xlim=c(-10, 10),ylim=c(0,1), main="")
for(iter in 1:nrow(theta.matrix)){
    ## record draw F ~ p(F | th,sig,y) (approx)
    f   <- fbar.H(xgrid,theta.matrix[iter,])
    lines(xgrid,f,col=iter,lty=3)
    fgrid <- rbind(fgrid,f)
  }
  ## add overall average (= posterior mean) to the plot
  fbar <- apply(fgrid,2,mean)
  lines(xgrid,fbar,lwd=3,col=2)
}

Posterior_G<- function(n, alpha, Gbase, theta){
  
}

Approximation 1st order

\[\theta_{n} \mid \theta_{n-1},..\theta_1 \sim \frac{k_n \alpha}{n} H + \frac{1}{n} \sum_i^{k_n} (n_i - \alpha)\delta_{\theta^*_j}\]

Approximation 2nd order

\[\theta_{n} \mid \theta_{n-1},..\theta_1 \sim \frac{k_n \alpha + \beta_n}{n + \beta_n} H + \frac{1}{n + \beta_n} \sum_i^{k_n} (n_i - \alpha)\delta_{\theta^*_j}\]

For PY process \(\beta_n = \theta\)

polya_urn_model <- function(base_measure, N_ball, alpha) {
  balls <- c()
  for (i in 1:N_ball) {
    if (runif(1) < alpha / (alpha + length(balls))) { 
      #Add a new ball color.
      new_color <- base_measure()
      balls <- c(balls, new_color)
        } else {
    #Pick out a ball from the urn, and add back a # ball of the same color
          
    ball <- balls[sample(1:length(balls), 1)] 
    balls <- c(balls, ball)
    } }
    balls 
}

alpha_vect <- c(1, 10, 100, 1000)
N_urns <- 3
sigma2 <- c(1,.4,.2)
N_draws <- 100
N_xaxis <- 200
x_axis <- seq(-3, 3, length = N_xaxis)
result <- NULL
  for (alpha in alpha_vect) {
    PU <- polya_urn_model(function() rnorm(1), N_draws, alpha)
    for (u in 1:N_urns) {
    res <- mapply(function(mean) dnorm(x_axis, rep(mean, N_xaxis),
                                       rep(sigma2[u], N_xaxis)), PU)
    res <- apply(res, 1, mean)
    new_draw <- cbind(res, x_axis, alpha, sigma2[u])
    result <- rbind(result, new_draw)
} }
result <- as.data.frame(result)
names(result) <- c("density", "x", "alpha", "sigma2")
DP_mixt <- qplot(data = result, y = density, x = x,
                 geom = c("line", "area")) +
  facet_grid(alpha ~ sigma2, labeller =
               label_bquote(rows = alpha == .(alpha),
                            cols = sigma^2 == .(sigma2))) +
  aes(color = as.factor(alpha)) +
  theme(legend.position = "none")
DP_mixt

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.